#include "Eigenstructure.h"

double Eigenstructure::getVectorInverse(int iK, int iP) const {
    return vectorsInverse[iP * K + iK];
}

double Eigenstructure::getPav(int ik, int iK, int m, int n) const {
    //pavs are stored in the raw data as matrix index-major order, go figure
    return pav[ik][(n + m*4) * K + iK];
}

double Eigenstructure::getx(int ik, int iK, int m, int n) const {
    return xRef(ik, iK, m, n);
}

double& Eigenstructure::xRef(int ik, int iK, int m, int n) const {
    return x[ik][n + m*4 + iK*16];
}

void Eigenstructure::setFromFile(int N, FILE* file) {
    this->N = N;
    K = 2*N+8;
    values = new double[K];
    vectorsInverse = new double[K*K];
    pav[0] = new double[K*16];
    pav[1] = new double[K*16];
    pav[2] = new double[K*16];

    fread(values, sizeof(double), K, file);
    fread(vectorsInverse, sizeof(double), K*K, file);
    fread(pav[0], sizeof(double), K*16, file);
    fread(pav[1], sizeof(double), K*16, file);
    fread(pav[2], sizeof(double), K*16, file);

    //convert from b-spline basis
    //based on http://lists.blender.org/pipermail/bf-blender-cvs/2007-October/011528.html
    //and http://www.codase.com/search/display?file=L2dlbnRvbzIvdmFyL3RtcC9yZXBvcy9jb2Rhc2UuYy9waXhpZS0xLjQuMS1yMS93b3JrL1BpeGllL3NyYy9wcmVjb21wL3ByZWNvbXAuY3Bw&lang=c%2B%2B
    x[0] = new double[K*16];
    x[1] = new double[K*16];
    x[2] = new double[K*16];

    double bsplineBasis[4][4] = {
        //N_0       N_1       N_2      N_3
        { 1.0/6.0,  4.0/6.0,  1.0/6.0, 0.0    }, // u^0
        {-3.0/6.0,  0.0    ,  3.0/6.0, 0.0    }, // u^1
        { 3.0/6.0, -6.0/6.0,  3.0/6.0, 0.0    }, // u^2
        {-1.0/6.0,  3.0/6.0, -3.0/6.0, 1.0/6.0}  // u^3
    };

    for (int ik = 0; ik < 3; ik++) {
        for (int iK = 0; iK < K; iK++) {
            //zero out the x matrix
            for (int xm = 0; xm < 4; xm++) for (int xn = 0; xn < 4; xn++) {
                xRef(ik, iK, xm, xn) = 0;
            }

            for (int pavm = 0; pavm < 4; pavm++) for (int pavn = 0; pavn < 4; pavn++) {
                double pavValue = getPav(ik, iK, pavm, pavn);
                //loop over indices of x (the elements of x being coefficients for u^xm v^xn)
                for (int xm = 0; xm < 4; xm++) for (int xn = 0; xn < 4; xn++)  {
                    xRef(ik, iK, xm, xn) += pavValue * bsplineBasis[xm][pavm] * bsplineBasis[xn][pavn];
                }
            }
        }
    }
}

Eigenstructure* Eigenstructure::readFile(int& Nmax) {
    FILE* file;
    Eigenstructure* result;
    file = fopen("ccdata50NT.dat", "rb");
    if (!file) return 0;

    fread ( &Nmax, sizeof(int), 1, file );
    result = new Eigenstructure[Nmax-2];

    for (int i = 0; i < Nmax-2 ; i++) {
        int N = i+3;
        result[i].setFromFile(N, file);
    }

    fclose(file);

    return result;
}

void Eigenstructure::print() {
    printf("N = %d\n", N);

    printf("eigenvalues:\n");

    for (int i = 0; i < K; i++) {
        printf("%6.3f ", values[i]);
    }
    printf("\n\n");
    printf("inverse eigenvectors:\n");

    for (int i = 0; i < K; i++) {
        for (int j = 0; j < K; j++) {
            printf("%6.3f ", vectorsInverse[i+j*K]);
        }
        printf("\n");
    }

    printf("\n");
    printf("pav:\n");

    for (int ik = 0; ik < 3; ik++) {
        printf("k=%d\n", ik);
        for (int iK = 0; iK < K; iK++) {
            for (int m = 0; m < 4; m++) {
                for (int n = 0; n < 4; n++) {
                    printf("%6.3f ", getPav(ik, iK, m, n));
                }
            }
            printf("\n");
        }
        printf("\n");
    }

    printf("x:\n");
    for (int ik = 0; ik < 3; ik++) {
        printf("k=%d\n", ik);
        for (int iK = 0; iK < K; iK++) {
            for (int m = 0; m < 4; m++) {
                for (int n = 0; n < 4; n++) {
                    printf("%6.3f ", xRef(ik, iK, m, n));
                }
            }
            printf("\n");
        }
        printf("\n");
    }
}

void Eigenstructure::printFile() {
    FILE* outfile = fopen("eigen.txt", "w");

    fprintf(outfile, "N = %d\n", N);

    fprintf(outfile, "eigenvalues:\n");

    for (int i = 0; i < K; i++) {
        fprintf(outfile, "%6.3f ", values[i]);
    }
    fprintf(outfile, "\n\n");
    fprintf(outfile, "inverse eigenvectors:\n");

    for (int i = 0; i < K; i++) {
        for (int j = 0; j < K; j++) {
            fprintf(outfile, "%6.3f ", vectorsInverse[i+j*K]);
        }
        fprintf(outfile, "\n");
    }

    fprintf(outfile, "\n");
    fprintf(outfile, "pav:\n");

    for (int ik = 0; ik < 3; ik++) {
        fprintf(outfile, "k=%d\n", ik);
        for (int iK = 0; iK < K; iK++) {
            for (int m = 0; m < 4; m++) {
                for (int n = 0; n < 4; n++) {
                    fprintf(outfile, "%6.3f ", getPav(ik, iK, m, n));
                }
            }
            fprintf(outfile, "\n");
        }
        fprintf(outfile, "\n");
    }

    fprintf(outfile, "x:\n");
    for (int ik = 0; ik < 3; ik++) {
        fprintf(outfile, "k=%d\n", ik);
        for (int iK = 0; iK < K; iK++) {
            for (int m = 0; m < 4; m++) {
                for (int n = 0; n < 4; n++) {
                    fprintf(outfile, "%6.3f ", xRef(ik, iK, m, n));
                }
            }
            fprintf(outfile, "\n");
        }
        fprintf(outfile, "\n");
    }

    fclose(outfile);
}
